Variation in vocal behavior under varying stress levels
Budgie INBRE stress and contact calls
Marcelo Araya-Salas, PhD & Timothy Wright, PhD
29-04-2022
Load packages
## add 'developer' to packages to be installed from github
x <- c("remotes", "lubridate", "readxl", "pbapply", "viridis", "ggplot2", "kableExtra", "knitr", "formatR", "MASS", "sp", "GGally", "brms", "lme4", "dplyr", "purrr", "forcats", "tidyr", "modelr", "tidybayes", "cowplot", "ggrepel", "posterior", "ggridges", "maRce10/PhenotypeSpace")
aa <- lapply(x, function(y) {
# get pakage name
pkg <- strsplit(y, "/")[[1]]
pkg <- pkg[length(pkg)]
# check if installed, if not then install
if (!pkg %in% installed.packages()[,"Package"]) {
if (grepl("/", y)) remotes::install_github(y, force = TRUE) else
install.packages(y)
}
# load package
a <- try(require(pkg, character.only = T), silent = T)
if (!a) remove.packages(pkg)
})Functions and global parameters
opts_knit$set(root.dir = "..")
# set evaluation false
opts_chunk$set(fig.width = 10, fig.height = 6, warning = FALSE, message = FALSE, tidy = TRUE)
read_excel_df <- function(...) data.frame(read_excel(...))
# for reading months in english format
sl <- Sys.setlocale(locale = "en_US.UTF-8")
standard_error <- function(x) sd(x)/sqrt(length(x))
cols <- viridis(10, alpha = 0.7)
# print results brms models
summary_brm_model <- function(x, gsub.pattern = NULL, gsub.replacement = NULL, xlab = "Effect size"){
summ <- summary(x)$fixed
fit <- x$fit
betas <- grep("^b_", names(fit@sim$samples[[1]]), value = TRUE)
hdis <- t(sapply(betas, function(y) hdi(fit@sim$samples[[1]][[y]]))
)
summ_table <- data.frame(summ, hdis, iterations = attr(fit, "stan_args")[[1]]$iter, chains = length(attr(fit, "stan_args")))
summ_table <- summ_table[rownames(summ_table) != "Intercept", c("Estimate", "Rhat", "Bulk_ESS", "l.95..CI", "u.95..CI", "iterations", "chains")]
summ_table <- as.data.frame(summ_table)
summ_table$Rhat <- round(summ_table$Rhat, digits = 3)
summ_table$CI_low <- round(unlist(summ_table$l.95..CI), digits = 3)
summ_table$CI_high <- round(unlist(summ_table$u.95..CI), digits = 3)
summ_table$l.95..CI <- summ_table$u.95..CI <- NULL
out <- lapply(betas, function(y) data.frame(variable = y, value = sort(fit@sim$samples[[1]][[y]], decreasing = FALSE)))
posteriors <- do.call(rbind, out)
posteriors <- posteriors[posteriors$variable != "b_Intercept", ]
if (!is.null(gsub.pattern) & !is.null(gsub.replacement))
posteriors$variable <- gsub(pattern = gsub.pattern, replacement = gsub.replacement, posteriors$variable)
gg <- ggplot(data = posteriors, aes(y = variable, x = value, fill = stat(quantile))) +
geom_density_ridges_gradient(quantile_lines = TRUE, quantile_fun = HDInterval::hdi, vline_linetype = 2) +
theme_classic() +
labs(x = xlab, y = "Predictor") +
xlim(range(c(min(summ_table[ , "CI_low"]), max(summ_table[ , "CI_high"])), 0)) +
geom_vline(xintercept = 0, col = "red") +
scale_fill_manual(values = c("transparent", "lightblue", "transparent"), guide = "none") + ggtitle(x$formula)
if (!is.null(gsub.pattern) & !is.null(gsub.replacement))
rownames(summ_table) <- gsub(pattern = gsub.pattern, replacement = gsub.replacement, rownames(summ_table))
summ_table$Rhat <- ifelse(summ_table$Rhat > 1.05, cell_spec(summ_table$Rhat, "html", color ="white", background = "red", bold = TRUE, font_size = 12), cell_spec(summ_table$Rhat, "html"))
signif <- summ_table[,"CI_low"] * summ_table[,"CI_high"] > 0
df1 <- kbl(summ_table, row.names = TRUE, escape = FALSE, format = "html", digits = 3)
df1 <- row_spec(kable_input = df1,row = which(summ_table$CI_low * summ_table$CI_high > 0), background = adjustcolor(cols[9], alpha.f = 0.3))
df1 <- kable_styling(df1, bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = FALSE, font_size = 12)
print(df1)
print(gg)
}# read ext sel tab calls
sels <- read.csv("./data/processed/tailored_budgie_calls_sel_tab.csv")
# keep only spectrographic parameters
sels <- sels[, c("sound.files", "selec", "duration", "meanfreq", "sd", "freq.median",
"freq.IQR", "time.IQR", "skew", "kurt", "sp.ent", "time.ent", "entropy", "meandom",
"mindom", "maxdom", "dfrange", "modindx", "meanpeakf")]
sels$ID <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][1])
sels$month <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][2])
sels$day <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][3])
sels$year <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][4])
sels$date <- paste(sels$day, substr(sels$month, 0, 3), sels$year, sep = "-")
sels$date <- as.Date(sels$date, format = "%d-%b-%Y")
# acoustic measurements
areas_by_week <- readRDS("./data/processed/acoustic_space_area_by_individual_and_week.RDS")
indiv_ovlp <- readRDS("./data/processed/acoustic_space_density_overlap_to_first_week_by_individual.RDS")
indiv_ovlp$treatment <- factor(indiv_ovlp$treatment, levels = c("Control", "Medium Stress",
"High Stress"))
group_ovlp <- readRDS("./data/processed/acoustic_space_density_overlap_to_group_by_week.RDS")
# group_ovlp <- do.call(rbind, lapply(unique(group_ovlp$ID), function(x) { X <-
# group_ovlp[group_ovlp$ID == x, ] X$overlap.to.group <- X$overlap.to.group -
# X$overlap.to.group[X$week == min(as.numeric(as.character(X$week)))]
# X$distance.to.group <- X$distance.to.group - X$distance.to.group[X$week ==
# min(as.numeric(as.character(X$week)))] return(X) }))
group_ovlp$treatment <- factor(group_ovlp$treatment, levels = c("Control", "Medium Stress",
"High Stress"))Counts per individual
df <- as.data.frame(table(sels$ID))
names(df) <- c("ID", "Sample_size")
df <- df[order(df$Sample_size, decreasing = FALSE), ]
kb <- kable(df, row.names = FALSE)
kb <- kable_styling(kb, bootstrap_options = c("striped", "hover", "condensed", "responsive"))
print(kb)| ID | Sample_size |
|---|---|
| 132YGMM | 6 |
| 125YGMM | 12 |
| 160YGHM | 16 |
| 310YGHM | 21 |
| 393YGHM | 25 |
| 303YGHM | 37 |
| 398WBHM | 41 |
| 365YGHM | 44 |
| 399YGLM | 46 |
| 300YGHM | 47 |
| 270YGHM | 64 |
| 407YGHM | 97 |
| 386WBHM | 100 |
| 377WWLM | 108 |
| 367WWNM | 119 |
| 371YYLM | 149 |
| 397YGHM | 175 |
| 378YYLM | 195 |
| 404WBHM | 196 |
| 258YGHM | 223 |
| 389WWLM | 230 |
| 262YGHM | 306 |
| 400YGHM | 360 |
| 388YGLM | 377 |
| 327YYLM | 404 |
| 396YBHM | 455 |
| 403WBHM | 566 |
| 356WBHM | 761 |
| 361YGHM | 767 |
| 402YGHM | 849 |
| 395WBHM | 854 |
| 360YGHM | 900 |
| 390WBHM | 935 |
| 405YBHM | 975 |
| 385YBHM | 1256 |
| 394YBHM | 1693 |
Add metadata
metadat <- read_excel_df("./data/raw/INBREStress_MasterDataSheet_14Nov19.xlsx")
# head(metadat)
sels$ID[sels$ID == "125YGMM"] <- "125YGHM"
sels$ID[sels$ID == "394YBHM"] <- "394WBHM"
# setdiff(sels$ID, metadat$Bird.ID) setdiff(metadat$Bird.ID, sels$ID)
sels$treatment <- sapply(1:nrow(sels), function(x) {
metadat$Treatment[metadat$Bird.ID == sels$ID[x]][1]
})
sels$treatment.room <- sapply(1:nrow(sels), function(x) {
metadat$Treatment.Room[metadat$Bird.ID == sels$ID[x]][1]
})
sels$round <- sapply(1:nrow(sels), function(x) {
metadat$Round[metadat$Bird.ID == sels$ID[x]][1]
})
sels$source.room <- sapply(1:nrow(sels), function(x) {
metadat$Source.Room[metadat$Bird.ID == sels$ID[x]][1]
})
sels$record.group <- sapply(1:nrow(sels), function(x) {
metadat$Record.Group[metadat$Bird.ID == sels$ID[x]][1]
})
# add week
out <- lapply(unique(sels$round), function(x) {
Y <- sels[sels$round == x, ]
min_date <- min(Y$date)
week_limits <- min_date + seq(0, 100, by = 7)
Y$week <- NA
for (i in 2:length(week_limits)) Y$week[Y$date >= week_limits[i - 1] & Y$date <
week_limits[i]] <- i - 1
return(Y)
})
sels <- do.call(rbind, out)
sels$cort.baseline <- sapply(1:nrow(sels), function(x) {
if (sels$week[x] == 1)
out <- metadat$D3.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 2)
out <- metadat$D7.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 3)
out <- metadat$D14.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 4)
out <- metadat$D21.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 5)
out <- metadat$D28.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]
return(out)
})
sels$cort.stress <- sapply(1:nrow(sels), function(x) {
if (sels$week[x] == 1)
out <- metadat$D3.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 2)
out <- metadat$D7.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 3)
out <- metadat$D14.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 4)
out <- metadat$D21.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 5)
out <- metadat$D28.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]
return(out)
})
sels$stress.response <- sels$cort.stress - sels$cort.baseline
sels$weight <- sapply(1:nrow(sels), function(x) {
if (sels$week[x] == 1)
out <- metadat$D3.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 2)
out <- metadat$D7.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 3)
out <- metadat$D14.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 4)
out <- metadat$D21.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 5)
out <- metadat$D28.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
return(out)
})
sels$breath.count <- sapply(1:nrow(sels), function(x) {
if (sels$week[x] == 1)
out <- metadat$D3.Breath.Count[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 2)
out <- metadat$D7.Breath.Count[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 3)
out <- metadat$D14.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 4)
out <- metadat$D21.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 5)
out <- metadat$D28.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
return(out)
})
# remove week 5
sels <- sels[sels$week != 5, ]agg_dat <- aggregate(selec ~ ID + week, data = sels, length)
# compare to week 1 agg_dat$call.count <- sapply(1:nrow(agg_dat), function(x) {
# baseline <- agg_dat$selec[agg_dat$week == 1 & agg_dat$ID == agg_dat$ID[x]] if
# (length(baseline) > 0) change <- agg_dat$selec[x] - baseline else change <-
# agg_dat$selec[x] return(change) } )
# without comparing to week 1
agg_dat$call.count <- sapply(1:nrow(agg_dat), function(x) agg_dat$selec[x])
agg_dat$selec <- NULL
agg_dat$baseline.CORT <- sapply(1:nrow(agg_dat), function(x) {
baseline <- sels$cort.baseline[sels$week == 1 & sels$ID == agg_dat$ID[x]]
current <- sels$cort.baseline[sels$week == agg_dat$week[x] & sels$ID == agg_dat$ID[x]]
if (length(baseline) > 0 & length(current) > 0)
change <- mean(current) - mean(baseline) else change <- NA
return(change)
})
agg_dat$stress.response <- sapply(1:nrow(agg_dat), function(x) {
baseline <- sels$stress.response[sels$week == 1 & sels$ID == agg_dat$ID[x]]
current <- sels$stress.response[sels$week == agg_dat$week[x] & sels$ID == agg_dat$ID[x]]
if (length(baseline) > 0 & length(current) > 0)
change <- mean(current) - mean(baseline) else change <- NA
return(change)
})
agg_dat$stress.CORT <- sapply(1:nrow(agg_dat), function(x) {
baseline <- sels$cort.stress[sels$week == 1 & sels$ID == agg_dat$ID[x]]
current <- sels$cort.stress[sels$week == agg_dat$week[x] & sels$ID == agg_dat$ID[x]]
if (length(baseline) > 0 & length(current) > 0)
change <- mean(current) - mean(baseline) else change <- NA
return(change)
})
agg_dat$weight <- sapply(1:nrow(agg_dat), function(x) {
baseline <- sels$weight[sels$week == 1 & sels$ID == agg_dat$ID[x]]
current <- sels$weight[sels$week == agg_dat$week[x] & sels$ID == agg_dat$ID[x]]
if (length(baseline) > 0 & length(current) > 0)
change <- mean(current) - mean(baseline) else change <- NA
return(change)
})
agg_dat$breath.rate <- sapply(1:nrow(agg_dat), function(x) {
baseline <- sels$breath.count[sels$week == 1 & sels$ID == agg_dat$ID[x]]
current <- sels$breath.count[sels$week == agg_dat$week[x] & sels$ID == agg_dat$ID[x]]
if (length(baseline) > 0 & length(current) > 0)
change <- mean(current) - mean(baseline) else change <- NA
return(change)
})
agg_dat$acoustic.area <- sapply(1:nrow(agg_dat), function(x) {
area <- areas_by_week$raref.area[areas_by_week$ID == agg_dat$ID[x] & areas_by_week$week ==
agg_dat$week[x]]
if (length(area) < 1)
area <- NA
return(area)
})
agg_dat$acoustic.distance <- sapply(1:nrow(agg_dat), function(x) {
distance <- indiv_ovlp$distance.to.first.week[indiv_ovlp$ID == agg_dat$ID[x] &
indiv_ovlp$week == agg_dat$week[x]]
if (length(distance) < 1)
distance <- NA
return(distance)
})
agg_dat$acoustic.overlap <- sapply(1:nrow(agg_dat), function(x) {
overlap <- indiv_ovlp$overlap.to.first.week[indiv_ovlp$ID == agg_dat$ID[x] &
indiv_ovlp$week == agg_dat$week[x]]
if (length(overlap) < 1)
overlap <- NA
return(overlap)
})
agg_dat$acoustic.overlap.to.group <- sapply(1:nrow(agg_dat), function(x) {
overlap <- group_ovlp$overlap.to.group[group_ovlp$ID == agg_dat$ID[x] & group_ovlp$week ==
agg_dat$week[x]]
if (length(overlap) < 1)
overlap <- NA
return(overlap)
})
agg_dat$treatment <- sapply(1:nrow(agg_dat), function(x) unique(sels$treatment[sels$ID ==
agg_dat$ID[x]]))
agg_dat$round <- sapply(1:nrow(agg_dat), function(x) unique(sels$round[sels$ID ==
agg_dat$ID[x]]))Physiological parameters
breath.count <- stack(metadat[, c("D3.Breath.Count", "D7.Breath.Count", "D14.Breath.Count",
"D21.Breath.Count", "D28.Breath.Count")])
weight <- stack(metadat[, c("D3.Bird.Weight..g.", "D7.Bird.Weight..g.", "D14.Bird.Weight..g.",
"D21.Bird.Weight..g.", "D28.Bird.Weight..g.")])
cort.stress <- stack(metadat[, c("D3.CORT.Stress", "D7.CORT.Stress", "D14.CORT.Stress",
"D21.CORT.Stress", "D28.CORT.Stress")])
cort.baseline <- stack(metadat[, c("D3.CORT.Baseline", "D7.CORT.Baseline", "D14.CORT.Baseline",
"D21.CORT.Baseline", "D28.CORT.Baseline")])
stress <- data.frame(metadat[, c("Bird.ID", "Treatment", "Round", "Treatment.Room")],
week = breath.count$ind, breath.count = breath.count$values, weight = weight$values,
cort.stress = cort.stress$values, cort.baseline = cort.baseline$values, stress.response = cort.stress$values -
cort.baseline$values)
# head(stress)
stress$week <- factor(sapply(strsplit(as.character(stress$week), "\\."), "[[", 1),
levels = c("D3", "D7", "D14", "D21", "D28"))
stress$day <- as.numeric(gsub("D", "", as.character(stress$week)))
stress$round <- paste("Round", stress$Round)
stress$treatment <- factor(stress$Treatment, levels = c("Control", "Medium Stress",
"High Stress"))
# remove 5th week
stress <- stress[stress$week != "D28", ]
stress_l <- lapply(stress$Bird.ID, function(x) {
X <- stress[stress$Bird.ID == x, ]
X$breath.count <- X$breath.count - X$breath.count[X$week == "D3"]
X$weight <- X$weight - X$weight[X$week == "D3"]
X$cort.stress <- X$cort.stress - X$cort.stress[X$week == "D3"]
X$cort.baseline <- X$cort.baseline - X$cort.baseline[X$week == "D3"]
X$stress.response <- X$stress.response - X$stress.response[X$week == "D3"]
return(X)
})
stress <- do.call(rbind, stress_l)
agg_stress <- aggregate(cbind(breath.count, weight, stress.response, cort.baseline) ~
week + treatment, stress, mean)
agg_stress_se <- aggregate(cbind(breath.count, weight, stress.response, cort.baseline) ~
week + treatment, stress, standard_error)
names(agg_stress_se) <- paste(names(agg_stress_se), ".se", sep = "")
agg_stress <- cbind(agg_stress, agg_stress_se[, 3:6])
agg_stress$Week <- 1:4
bs <- 22
gg_breath.count <- ggplot(data = agg_stress, aes(x = Week, y = breath.count, fill = treatment)) +
geom_bar(stat = "identity") + geom_errorbar(aes(ymin = breath.count - breath.count.se,
ymax = breath.count + breath.count.se), width = 0.1) + scale_fill_viridis_d(begin = 0.1,
end = 0.9) + facet_wrap(~treatment, ncol = 3, scale = "fixed") + labs(y = "Mean change in breath\nrate (breaths/min)",
x = "Week") + theme_classic(base_size = bs) + theme(legend.position = "none")
gg_weight <- ggplot(data = agg_stress, aes(x = Week, y = weight, fill = treatment)) +
geom_bar(stat = "identity") + geom_errorbar(aes(ymin = weight - weight.se, ymax = weight +
weight.se), width = 0.1) + scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~treatment,
ncol = 3, scale = "fixed") + labs(y = "Mean change in \nweight (grams)", x = "Week") +
theme_classic(base_size = bs) + theme(legend.position = "none")
gg_cort.baseline <- ggplot(data = agg_stress, aes(x = Week, y = cort.baseline, fill = treatment)) +
geom_bar(stat = "identity") + geom_errorbar(aes(ymin = cort.baseline - cort.baseline.se,
ymax = cort.baseline + cort.baseline.se), width = 0.1) + scale_fill_viridis_d(begin = 0.1,
end = 0.9) + facet_wrap(~treatment, ncol = 3, scale = "fixed") + labs(y = "Mean change in\nbaseline CORT (ng/mL)",
x = "Week") + theme_classic(base_size = bs) + theme(legend.position = "none")
gg_stress.response <- ggplot(data = agg_stress, aes(x = Week, y = stress.response,
fill = treatment)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = stress.response -
stress.response.se, ymax = stress.response + stress.response.se), width = 0.1) +
scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~treatment, ncol = 3,
scale = "fixed") + labs(y = "Mean change in stress\nresponse CORT (ng/mL)", x = "Week") +
theme_classic(base_size = bs) + theme(legend.position = "none")
cowplot::plot_grid(gg_breath.count, gg_weight, gg_cort.baseline, gg_stress.response)
Stats
Models: Predicted physio measure ~ treatment + week (continuous) + IndRandom
Variables (Difference from Week 1): weight, BR, baseline CORT, Stress CORT, Stress Response
responses <- c("baseline.CORT", "stress.response", "stress.CORT", "weight", "breath.rate")
predictors <- c("~ treatment + week + (1|ID) + (1|round)")
formulas <- expand.grid(responses = responses, predictors = predictors, stringsAsFactors = FALSE)
vars_to_scale <- c(responses, "week")
# remove week 1
sub_agg_dat <- agg_dat[agg_dat$week != 1, ]
for (i in vars_to_scale) sub_agg_dat[, vars_to_scale] <- scale(sub_agg_dat[, vars_to_scale])
physio_models <- lapply(1:nrow(formulas), function(x) {
sub_dat <- sub_agg_dat[!is.na(sub_agg_dat[names(sub_agg_dat) == formulas$responses[x]]),
]
sub_dat
mod <- try(brm(formula = paste(formulas$responses[x], formulas$predictors[x]),
iter = 20000, silent = 2, data = sub_dat, control = list(adapt_delta = 0.9),
chains = 4), silent = TRUE)
return(mod)
})
names(physio_models) <- paste(formulas$responses, formulas$predictors)
saveRDS(physio_models, "./data/processed/physiological_response_models.RDS")physio_models <- readRDS("./data/processed/physiological_response_models.RDS")
ggl <- list()
for (x in 1:length(physio_models)) {
cat(paste("- ", names(physio_models)[x], "\n"))
ggl[[length(ggl) + 1]] <- summary_brm_model(physio_models[[x]], gsub.pattern = "b_treatment|b_",
gsub.replacement = "")
}- baseline.CORT ~ treatment + week + (1|ID) + (1|round)
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high treatmentHighStress 0.763 1.001 10563.66 20000 4 0.114 1.417 treatmentMediumStress 0.180 1 11354.05 20000 4 -0.511 0.868 week 0.000 1 30967.77 20000 4 -0.142 0.142 
- stress.response ~ treatment + week + (1|ID) + (1|round)
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high treatmentHighStress -0.493 1.001 8754.839 20000 4 -1.182 0.191 treatmentMediumStress 0.002 1.001 7335.265 20000 4 -0.768 0.760 week -0.123 1.001 38086.430 20000 4 -0.255 0.007 
- stress.CORT ~ treatment + week + (1|ID) + (1|round)
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high treatmentHighStress -0.217 1 12610.60 20000 4 -0.908 0.469 treatmentMediumStress 0.076 1 13662.12 20000 4 -0.697 0.837 week -0.153 1 37967.80 20000 4 -0.285 -0.020 
- weight ~ treatment + week + (1|ID) + (1|round)
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high treatmentHighStress -0.315 1 10029.64 20000 4 -0.966 0.343 treatmentMediumStress -0.126 1 10793.59 20000 4 -0.835 0.587 week -0.344 1 21762.91 20000 4 -0.478 -0.210 
- breath.rate ~ treatment + week + (1|ID) + (1|round)
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high treatmentHighStress 0.145 1.001 21047.468 20000 4 -0.158 0.438 treatmentMediumStress 0.058 1.001 8702.607 20000 4 -0.275 0.380 week -0.785 1 35693.191 20000 4 -0.898 -0.670 
names(ggl) <- names(physio_models)
Barplot and effect sizes side-by-side
cowplot::plot_grid(gg_breath.count + theme_classic(base_size = 10) + theme(legend.position = "none"),
ggl[[grep("breath", names(ggl))]])
cowplot::plot_grid(gg_weight + theme_classic(base_size = 10) + theme(legend.position = "none"),
ggl[[grep("weight", names(ggl))]])
cowplot::plot_grid(gg_cort.baseline + theme_classic(base_size = 10) + theme(legend.position = "none"),
ggl[[grep("base", names(ggl))]])
cowplot::plot_grid(gg_stress.response + theme_classic(base_size = 10) + theme(legend.position = "none"),
ggl[[grep("response", names(ggl))]])
Takeaways
Breath rate decreases gradually with time across after the first week
Stress response is higher in “high stress” birds compared to first week
Acoustic space projection
t-SNE
scale_param <- scale(sels[, c("duration", "meanfreq", "sd", "freq.median", "freq.IQR",
"time.IQR", "skew", "kurt", "sp.ent", "time.ent", "entropy", "meandom", "mindom",
"maxdom", "dfrange", "modindx", "meanpeakf")])
tsne <- Rtsne(scale_param, dims = 2, perplexity = 30, verbose = FALSE, max_iter = 5000)
saveRDS(tsne, "./data/processed/tsne_on_acoustic_parameters_jun_2021.RDS")tsne <- readRDS("./data/processed/tsne_on_acoustic_parameters_jun_2021.RDS")
Y <- as.data.frame(tsne$Y)
names(Y) <- c("TSNE1", "TSNE2")
sels <- data.frame(sels, Y)
sels$treatment <- factor(sels$treatment, levels = c("Control", "Medium Stress", "High Stress"))
ggplot(sels, aes(x = TSNE1, y = TSNE2, col = as.factor(treatment))) + geom_point() +
labs(color = "Treatment") + scale_color_viridis_d(alpha = 0.4) + theme_classic(base_size = 25) +
guides(colour = guide_legend(override.aes = list(size = 10)))
Behavioral parameters
agg_call.count <- aggregate(cbind(call.count, acoustic.overlap.to.group) ~ week +
treatment, agg_dat, mean)
agg_behav <- aggregate(cbind(acoustic.area, acoustic.overlap) ~ week + treatment,
agg_dat, mean)
agg_call.count_se <- aggregate(cbind(call.count, acoustic.overlap.to.group) ~ week +
treatment, agg_dat, standard_error)
agg_behav_se <- aggregate(cbind(acoustic.area, acoustic.overlap) ~ week + treatment,
agg_dat, standard_error)
agg_behav_se <- merge(agg_call.count_se, agg_behav_se, all = TRUE)
names(agg_behav_se) <- paste(names(agg_behav_se), ".se", sep = "")
agg_behav <- merge(agg_call.count, agg_behav, all = TRUE)
agg_behav <- cbind(agg_behav, agg_behav_se[, 3:6])
bs <- 22
agg_behav$treatment <- factor(agg_behav$treatment, levels = c("Control", "Medium Stress",
"High Stress"))
gg_call.count <- ggplot(data = agg_behav, aes(x = week, y = call.count, fill = treatment)) +
geom_bar(stat = "identity") + geom_errorbar(aes(ymin = call.count - call.count.se,
ymax = call.count + call.count.se), width = 0.1) + scale_fill_viridis_d(begin = 0.1,
end = 0.9) + facet_wrap(~treatment, ncol = 3, scale = "fixed") + labs(y = "Vocal output (number of calls)",
x = "Week") + theme_classic(base_size = bs) + theme(legend.position = "none")
gg_acoustic.area <- ggplot(data = agg_behav, aes(x = week, y = acoustic.area, fill = treatment)) +
geom_bar(stat = "identity") + geom_errorbar(aes(ymin = acoustic.area - acoustic.area.se,
ymax = acoustic.area + acoustic.area.se), width = 0.1) + scale_fill_viridis_d(begin = 0.1,
end = 0.9) + facet_wrap(~treatment, ncol = 3, scale = "fixed") + labs(y = "Mean change in \n acoustic area",
x = "Week") + theme_classic(base_size = bs) + theme(legend.position = "none")
gg_acoustic.overlap <- ggplot(data = agg_behav, aes(x = week, y = acoustic.overlap,
fill = treatment)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = acoustic.overlap -
acoustic.overlap.se, ymax = acoustic.overlap + acoustic.overlap.se), width = 0.1) +
scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~treatment, ncol = 3,
scale = "fixed") + labs(y = "Acoustic overlap to self in week 1", x = "Week") +
theme_classic(base_size = bs) + theme(legend.position = "none")
gg_acoustic.overlap.group <- ggplot(data = agg_behav, aes(x = week, y = acoustic.overlap.to.group,
fill = treatment)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = acoustic.overlap.to.group -
acoustic.overlap.to.group.se, ymax = acoustic.overlap.to.group + acoustic.overlap.to.group.se),
width = 0.1) + scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~treatment,
ncol = 3, scale = "fixed") + labs(y = "Acoustic overlap to group", x = "Week") +
theme_classic(base_size = bs) + theme(legend.position = "none")
cowplot::plot_grid(gg_call.count, gg_acoustic.area, gg_acoustic.overlap, gg_acoustic.overlap.group)
Stats
Model: Predicted behavior ~ treatment + week (continuous) + IndRandom
Variables: # calls, Distance moved from self in first week, Overlap to original acoustic space, Match to group repertoire, Maybe overall size of acoustic space
Do as comparison to week one using rarefacted calls and kernel density
responses <- c("call.count", "acoustic.area", "acoustic.overlap", "acoustic.overlap.to.group")
predictors <- c("~ treatment + week + (1|ID) + (1|round)")
formulas <- expand.grid(responses = responses, predictors = predictors, stringsAsFactors = FALSE)
vars_to_scale <- c(responses, "week")
for (i in vars_to_scale) agg_dat[, vars_to_scale] <- scale(agg_dat[, vars_to_scale])
behav_models <- lapply(1:nrow(formulas), function(x) {
sub_dat <- agg_dat[!is.na(agg_dat[names(agg_dat) == formulas$responses[x]]),
]
# remove week 1
if (!grepl("count|group", formulas$responses[x]))
sub_dat <- sub_dat[sub_dat$week != 1, ]
mod <- try(brm(formula = paste(formulas$responses[x], formulas$predictors[x]),
iter = 50000, silent = 2, data = sub_dat, control = list(adapt_delta = 0.9),
chains = 4), silent = TRUE)
return(mod)
})
names(behav_models) <- paste(formulas$responses, formulas$predictors)
saveRDS(behav_models, "./data/processed/behavioral_response_models.RDS")behav_models <- readRDS("./data/processed/behavioral_response_models.RDS")
ggl <- list()
for (x in 1:length(behav_models)) {
cat(paste("- ", names(behav_models)[x], "\n \n---\n "))
ggl[[length(ggl) + 1]] <- summary_brm_model(behav_models[[x]], gsub.pattern = "b_treatment|b_",
gsub.replacement = "")
}- call.count ~ treatment + week + (1|ID) + (1|round) —
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high treatmentHighStress 0.790 1.001 11548.09 50000 4 0.213 1.380 treatmentMediumStress 0.398 1.001 17599.85 50000 4 -0.182 0.985 week -0.059 1 102691.50 50000 4 -0.203 0.086 
- acoustic.area ~ treatment + week + (1|ID) + (1|round) —
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high treatmentHighStress -0.206 1.001 10719.216 50000 4 -1.059 0.624 treatmentMediumStress -0.708 1 18329.938 50000 4 -1.567 0.143 week 0.035 1.003 1444.876 50000 4 -0.094 0.152 
- acoustic.overlap ~ treatment + week + (1|ID) + (1|round) —
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high treatmentHighStress 0.910 1 25225.85 50000 4 0.065 1.757 treatmentMediumStress 0.233 1 25160.95 50000 4 -0.699 1.151 week -0.169 1 76169.60 50000 4 -0.307 -0.032 
- acoustic.overlap.to.group ~ treatment + week + (1|ID) + (1|round) —
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high treatmentHighStress 0.351 1.001 28237.431 50000 4 -0.292 1.015 treatmentMediumStress 0.380 1.002 3209.631 50000 4 -0.280 1.046 week -0.223 1.001 41111.148 50000 4 -0.385 -0.057 
names(ggl) <- names(behav_models)
Barplot and effect sizes side-by-side
cowplot::plot_grid(gg_call.count + theme_classic(base_size = 10) + theme(legend.position = "none"),
ggl[[grep("count", names(ggl))]])
cowplot::plot_grid(gg_acoustic.area + theme_classic(base_size = 10) + theme(legend.position = "none"),
ggl[[grep("area", names(ggl))]])
cowplot::plot_grid(gg_acoustic.overlap + theme_classic(base_size = 10) + theme(legend.position = "none"),
ggl[[grep("overlap ~", names(ggl))]])
cowplot::plot_grid(gg_acoustic.overlap.group + theme_classic(base_size = 10) + theme(legend.position = "none"),
ggl[[grep("group", names(ggl))]])
Takeaways
Higher vocal output in “high stress” birds compared to control
Higher acoustic overlap to themselves in week 1 for “high stress” birds compared to control
Decrease in overlap to themselves in week 1 across time
Session information
## R version 4.1.1 (2021-08-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=es_ES.UTF-8
## [7] LC_PAPER=es_CR.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] PhenotypeSpace_0.1.0 warbleR_1.1.27 NatureSounds_1.0.4
## [4] seewave_2.2.0 tuneR_1.4.0 ggridges_0.5.3
## [7] posterior_1.2.1 ggrepel_0.9.1 cowplot_1.1.1
## [10] tidybayes_3.0.2 modelr_0.1.8 tidyr_1.2.0
## [13] forcats_0.5.1 purrr_0.3.4 dplyr_1.0.8
## [16] lme4_1.1-29 Matrix_1.3-4 brms_2.17.0
## [19] Rcpp_1.0.8.3 GGally_2.1.2 sp_1.4-6
## [22] MASS_7.3-54 formatR_1.12 knitr_1.38
## [25] kableExtra_1.3.4 ggplot2_3.3.5 viridis_0.6.2
## [28] viridisLite_0.4.0 pbapply_1.5-0 readxl_1.4.0
## [31] lubridate_1.8.0 remotes_2.4.2
##
## loaded via a namespace (and not attached):
## [1] utf8_1.2.2 tidyselect_1.1.2 fftw_1.0-7
## [4] htmlwidgets_1.5.4 grid_4.1.1 munsell_0.5.0
## [7] codetools_0.2-18 DT_0.22 miniUI_0.1.1.1
## [10] withr_2.5.0 Brobdingnag_1.2-7 spatstat.random_2.2-0
## [13] colorspace_2.0-3 highr_0.9 rstudioapi_0.13
## [16] stats4_4.1.1 dtw_1.22-3 tensor_1.5
## [19] bayesplot_1.9.0 labeling_0.4.2 rstan_2.21.5
## [22] polyclip_1.10-0 farver_2.1.0 bridgesampling_1.1-2
## [25] coda_0.19-4 vctrs_0.4.1 generics_0.1.2
## [28] xfun_0.30 R6_2.5.1 markdown_1.1
## [31] HDInterval_0.2.2 gamm4_0.2-6 projpred_2.0.2
## [34] bitops_1.0-7 spatstat.utils_2.3-0 reshape_0.8.9
## [37] assertthat_0.2.1 promises_1.2.0.1 scales_1.2.0
## [40] gtable_0.3.0 processx_3.5.3 goftest_1.2-3
## [43] rlang_1.0.2 systemfonts_1.0.4 splines_4.1.1
## [46] spatstat.geom_2.4-0 broom_0.8.0 checkmate_2.0.0
## [49] inline_0.3.19 yaml_2.3.5 reshape2_1.4.4
## [52] abind_1.4-5 threejs_0.3.3 crosstalk_1.2.0
## [55] backports_1.4.1 httpuv_1.6.5 tensorA_0.36.2
## [58] tools_4.1.1 ellipsis_0.3.2 spatstat.core_2.4-2
## [61] raster_3.5-15 jquerylib_0.1.4 RColorBrewer_1.1-3
## [64] proxy_0.4-26 plyr_1.8.7 base64enc_0.1-3
## [67] RCurl_1.98-1.6 ps_1.6.0 prettyunits_1.1.1
## [70] rpart_4.1-15 deldir_1.0-6 zoo_1.8-10
## [73] cluster_2.1.2 magrittr_2.0.3 ggdist_3.1.1
## [76] colourpicker_1.1.1 mvtnorm_1.1-3 matrixStats_0.62.0
## [79] shinyjs_2.1.0 mime_0.12 evaluate_0.15
## [82] arrayhelpers_1.1-0 xtable_1.8-4 shinystan_2.6.0
## [85] gridExtra_2.3 rstantools_2.2.0 compiler_4.1.1
## [88] tibble_3.1.6 crayon_1.5.1 minqa_1.2.4
## [91] StanHeaders_2.21.0-7 htmltools_0.5.2 mgcv_1.8-36
## [94] later_1.3.0 RcppParallel_5.1.5 DBI_1.1.1
## [97] boot_1.3-28 permute_0.9-7 cli_3.2.0
## [100] parallel_4.1.1 igraph_1.3.0 pkgconfig_2.0.3
## [103] signal_0.7-7 terra_1.5-21 spatstat.sparse_2.1-1
## [106] xml2_1.3.3 svUnit_1.0.6 dygraphs_1.1.1.6
## [109] svglite_2.1.0 bslib_0.3.1 webshot_0.5.3
## [112] rvest_1.0.2 stringr_1.4.0 distributional_0.3.0
## [115] callr_3.7.0 digest_0.6.29 vegan_2.6-2
## [118] spatstat.data_2.2-0 rmarkdown_2.13 cellranger_1.1.0
## [121] shiny_1.7.1 gtools_3.9.2 rjson_0.2.21
## [124] nloptr_2.0.0 lifecycle_1.0.1 nlme_3.1-152
## [127] jsonlite_1.8.0 fansi_1.0.3 pillar_1.7.0
## [130] lattice_0.20-44 loo_2.5.1 fastmap_1.1.0
## [133] httr_1.4.2 pkgbuild_1.3.1 glue_1.6.2
## [136] xts_0.12.1 shinythemes_1.2.0 stringi_1.7.6
## [139] sass_0.4.1